Graphite— diamond phase coexistence study employing a neural-network mapping of 

the ab initio potential energy surface 



O 

(N 

S-H . 

< 

O 

(N 



Rustam Z. Khaliullin/'B Hagai Eshet, 1 Thomas D. Kiihne, 1,2 Jorg Behler, 3 and Michele Parrinello 1 

Department of Chemistry and Applied Biosciences, ETH Zurich, 
USI Campus, via G. Buffi 13, 6900 Lugano, Switzerland 
2 Department of Physics and Division of Engineering and Applied Sciences, 
Harvard University, Cambridge, Massachusetts 02138, USA 
3 Lehrstuhl fur Theoretische Chemie, Ruhr-Universitat Bochum, D-44780 Bochum, Germany 

(Dated: April 21, 2010) 

An interatomic potential for the diamond and graphite phases of carbon has been created using 
a neural-network (NN) representation of the ab initio potential energy surface. The NN potential 
combines the accuracy of a first-principle description of both phases with the efficiency of empirical 
force fields and allows one to perform, for the first time, a molecular dynamics study, of ab initio 
quality, of the thermodynamics of graphite-diamond coexistence. Good agreement between the 
experimental and calculated coexistence curves is achieved if nuclear quantum effects are included 
in the simulation. 
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The ability of carbon atoms to form strong chemical 
bonds with a variety of coordination numbers leads to a 
remarkably wide range of physical properties of the con- 
densed phases of carbon. The diamond phase is a three- 
dimensional network of four-fold coordinated atoms char- 
acterized by a very low electrical conductivity and ex- 
treme hardness. Unlike diamond, the graphite phase is 
semimetallic and made up of planes of three-fold coordi- 
nated atoms. It behaves as a lubricant because of weak 
van der Waals (vdW) bonding between the planes. 

In spite of the great fundamental and practical im- 
portance of graphite and diamond the characterization 
of these phases and their mutual transformation is far 
from complete especially in the region of high pressures 
and temperatures which are difficult to access experimen- 
tally. Although computer simulations based on density 
functional theory (DFT) provide a comprehensive frame- 
work for modeling a variety of carbon polymorphs, they 
become computationally too demanding for the gener- 
ation of long molecular dynamics (MD) trajectories for 
large systems (nanosecond-long trajectories are required 
to study thermodynamics and mechanism of phases tran- 
sitions). On the other hand, the construction of ac- 
curate and computationally efficient potentials capable 
of describing the wide range of interactions in carbon 
is still a challenge. Many simple force fields devel- 
oped for covalent systems such as the embedded atom 
method the Stillinger- Weber approach Q, and the 
bond-order potential of Tersoff @ have only limited suc- 
cess in modeling carbon phases. More elaborate poten- 
tials such as the Brenner potential the environment- 
dependent interaction potential [j| , and a family of long- 
range carbon bond-order potentials @ significantly im- 
prove the description of carbon structures by incorporat- 
ing 7r-bonding effects and vdW interactions. Neverthe- 
less, even the most sophisticated empirical potentials do 



not always give a correct description of all properties or 
phenomena of interest. 

In the present paper, we followed a different approach 
for modeling solid phases of carbon such as diamond and 
graphite. Instead of representing the interatomic interac- 
tion energy by an analytic function fitted to experimen- 
tal (or ab initio) data we created an accurate mapping of 
the relevant portion of the ab initio potential energy sur- 
face (PES) using a recently developed high-dimensional 
neural network (NN) approach Q • This approach elimi- 
nates the requirement to guess a complicated functional 
form for the interatomic potential. Accurate mapping 
ensures that all properties determined by the topology 
of the PES are described with an accuracy comparable 
with that of DFT. Furthermore, PES mapping allows 
one to examine nuclear quantum effects in MD simula- 
tions from first principles whereas empirical potentials 
attempt to incorporate such effects through parameter- 
ization. From a computational standpoint, the NN en- 
ergies, forces, and stress tensor are evaluated with the 
speed of empirical potentials [g, l9j thus enabling an MD 
study of graphite-diamond coexistence of unprecedented 
accuracy. 

Neural networks have been successfully used to inter- 
polate the PES of simple chemical systems for the last 
decade 

However, an NN-based method that can 
be used to map the high-dimensional PES of bulk systems 
and large clusters has been introduced only recently @- 
9]. This mapping of the ab initio PES is performed by 
optimizing NN parameters to reproduce the ab initio en- 
ergies of many thousands of structures in a training set. 
The overfitting (i.e. obtaining a good fit to the training 
data, but performing less accurately when making predic- 
tions) is controlled by testing the performance of the NN 
for an independent test set not used in the optimization. 

The accuracy of the reference ab initio energies is of 
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paramount importance while training the network. It 
is known that conventional local- and semilocal den- 
sity functionals cannot describe the long-range electron 
correlations that are responsible for the vdW interac- 
tions between graphite sheets [l^]. To account for the 
dispersion forces in graphite, we employed the Perdew- 
Burke-Ernzerhof (PBE) functional in combination with 
the dispersion corrected atom centered pseudopotential 
(DCACP) 16], which has been shown to perform well 
for graphene sheets and aromatic compounds [3, f]~7j j . 
Extensive tests were performed to demonstrate that 
DCACP closely reproduces the experimental lattice con- 
stants as well as elastic and vibrational properties of di- 
amond and graphite (TABLE I}. The ABINIT pack- 
age [Hj was used to perform the ab initio calculations. 
A dense mesh of fc-points and a large plane-wave cutoff 
of 170 Ry were used for all structures so as to ensure 
convergence of the total energy to 1 meV/atom. 

The initial fitting of the carbon NN potential was 
performed on crystal structures that included the zero- 
temperature and randomly distorted structures of cubic 
and hexagonal diamond, hexagonal and rhombohedral 
graphite in the pressure range from -10 to 200 GPa. 
After the initial training, the NN was improved self- 
consistently by iterative repetition of the NN-driven 
MD simulations, collection of new structures emerging 
from the simulations, calculation of the DFT energies 
for the physically relevant structures, and refinement 
of the NN. These iterations were performed until the 
root mean squared error (RMSE) of the new structures 
not included in the fit converged to the RMSE of the 
test set. After the self-consistent procedure the DFT 
dataset contained ~60,000 DFT energies corresponding 
to more than 700,000 atomic environments. 10% of all 
structures were randomly chosen for the test set. The 
best fit was obtained for a NN with 2 hidden layers, 
each of which contains 25 nodes (the total number of 
the NN parameters is 1901). The RMSE of the train- 
ing set is 4.0 meV/atom, while the RMSE of the test 
set is 4.9 meV/atom. The maximum absolute errors 
are 41.5 meV/atom and 46.7 meV/atom for the train- 
ing and test sets, respectively. The largest errors are 
attributed to highly distorted graphite structures that 
are accessible only at temperatures of 4000-5000 K. At 
these high temperatures the errors are small compared 
to ksT ~ 340 — 430 meV so the quality of the relevant 
ensemble averages is essentially maintained. 

To check the accuracy of the NN potential we calcu- 
lated lattice constants, stiffness coefficients, and vibra- 
tional frequencies for the zero-temperature structures of 
cubic diamond and hexagonal graphite. The lattice con- 
stants were determined by minimization of the NN po- 
tential energy fitted using the Murnaghan equation [26] 
in the case of diamond and by a two-dimensional fourth 
order polynomial in the case of graphite. The second- 
order elastic constants were calculated by fitting the en- 
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FIG. 1: NN prediction and experimental data [29|] for the 
Co lattice parameter of hexagonal graphite as a function of 
pressure. T = 295 K. The NN values are computed from 
constant-pressure MD simulations with a quantum Langevin 
thermostat 13011 - 



ergy as a function of an appropriate cell distortion to 
a parabola (27L l28| while vibrational frequencies were 
obtained by diagonalizing the dynamical matrix. The 
computed quantities are summarized and compared with 
DFT and experimental values in TABLElB The NN accu- 
rately reproduces DFT results for structural, elastic and 
vibrational properties of diamond [IH . All properties of 
graphite determined by the strong in-plane interactions 
(ao, en, C12, Tzo) are also described well by the NN. 
However, the relative error between DFT and NN values 
is generally larger for the properties determined by weak 
interplanar interactions (e.g. C33, C44, C13). Nevertheless, 
the NN description of one of the most important struc- 
tural characteristics of graphite - the interlayer distance 
- is remarkably accurate for a wide range of pressures 
(FIG. Q]). 

The graphite-diamond coexistence line was determined 
by locating points of equal chemical potential in the P-T 
plane. This was done in three steps. First, we calculated 
the Helmholz free energy F NN (T , p ) of both phases at 
To = 2000 K by thermodynamic integration using Ein- 
stein crystals as the reference systems [3l[ 

F NN (T ,po) = F™(T , P0 ) + J\^),d\, (1) 

where U x = XU NN + (1 - X)U EIN . 

In the next step, the chemical potentials were evalu- 
ated by integrating the free energy as a function of den- 
sity starting from po 32] 



M NN (T ,p) = lf NN (T , Po ) + ^ + 6(T )ln^ 
N po po 

+ b(To) + c(To){2p - po). 



(2) 



Parameters a(Tn), b(To), and c(Tq) were determined by 
fitting the pressure dependence on density using 



P(T , p) = a(T ) + b{To)p + c(T )p 2 



(3) 
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TABLE I: Structural, elastic, and vibrational properties of graphite and diamond. 



Lattice const. (A) Elastic constants (GPa) Freq. (cm 1 ) 



Hex. graphite 


ao 


Co 


Bo 


Cll 


Cl2 


C33 


C44 Cl3 


r.zo 


r LO /TO 


PBE a 


2.461 


8.712 


2.4 


1240 


a 


2.4 


-0.5 


1561, 1561 


881 


PBE, DCACP 


2.467 


6.815 


37 


1069 


162 


40 


5 -4 


1553, 1573 


870 


NN 


2.467 


6.688 


48 


1080 


179 


52 


7 


1527, 1530 


834 


Exp. 6 


2.461 


6.705 


36.4±1.1 


1060±16 


180±20 


36.5±1 


4.0±0.4 15±5 


1575 


861 


Cub. diamond 


ao 




B 


Cll 


Cl2 


c 44 




To 




PBE a 


3.568 




432 


1060 


125 


562 




1289 




PBE, DCACP 


3.570 




439 


1056 


130 


567 




1292 




NN 


3.569 




434 


1016 


142 


580 




1295 




Exp. c 


3.567 




442 


1076.4±0.2 


125.2±2.3 


577.4±1.4 




1332 





"Results of calculations with the standard Vanderbilt ultrasoft PP from Ref. Il5l cn + C12 value from Ref. I IE 



'Lattice constants from Ref. 
^Lattice constants from Ref. 



clastic constants from Ref. 
clastic constants from Ref. 



vibrational frequencies from Refsj2 
vibrational frequency from Ref. |2 



Finally, the coexistence line was traced by integrating 
the Clausius-Clapeyron equation ^ = ^f{AV) usm S the 
predictor-corrector scheme of Kofke [33j . 

It is important to emphasize that long MD trajecto- 
ries are essential to obtain statistically accurate results 
for all three steps. Furthermore, it is desirable to per- 
form simulations using large systems as finite size effects 
can introduce significant errors to the free energies eval- 
uated by thermodynamic integration [34l ]. Hence, direct 
ab initio MD simulations for large systems (especially 
with a large plane wave cutoff and a dense fc-point mesh) 
are computationally very demanding for the evaluation 
of free energies, whereas, the NN provides an affordable 
and accurate method to determine the coexistence line. 

NN-driven MD simulations were performed for 512 
atoms of cubic diamond (cubic box, p = 173.94 ran -3 ) 
and 960 atoms of hexagonal graphite (4 layers, cell size 
ratio 2.024:2.104:1, p = 120.02 nuT 3 ). The temperature 
was controlled using a colored-noise Langevin thermostat 
that was tuned to provide the optimum sampling effi- 
ciency over all relevant vibrational modes 35] . The time 
step was set to 0.7 fs. The integral in Eq. Q] was evalu- 
ated numerically by the Gauss-Legendre quadrature with 
20 points. At each value of A, the average value of the 
integrand and its statistical error were obtained from a 
133 ps trajectory. State points along the 2000 K isotherm 
were obtained from NPT simulations governed by Nose- 
Hoover equations of motion with Langevin noise on the 
particle and cell velocities [35l . [3(3]. Averaging over a 
95 ps trajectory was performed for each state point. The 
predictor-corrector algorithm was iterated until pressure 
had converged to less than 0.05 GPa that required 2-3 it- 
erations of 50 ps each. The total simulation time required 
to obtain the coexistence line totals ~5 ns for each phase 
clearly demonstrating the advantage of the NN approach 
in comparison with the direct ab initio simulation. 

We performed two separate calculations of the coexis- 
tence line. In the first simulation, the Langevin thermo- 
stat was tuned to reproduce quantum-mechanical behav- 



ior of carbon nuclei using a recently published method of 
Ceriotti et al. (30j In the second simulation, the thermo- 
stat was fitted to obtain classical behavior of the nuclei. 

Two graphite-diamond coexistence lines determined as 
the intersection of the /x NN (T, P(p)) planes in classical 
and quantum simulations are shown in FIG. [2] We veri- 
fied that the coexistence lines are calculated correctly by 
independent thermodynamic integration at To = 300 K 
and T = 1000 K (indicated by red points in FIG. 0}. 
Comparison with the experimental data 37, [38[ in the 
temperature interval from 1500 to 3000 K reveals that 
the NN overestimates the transition pressure by approxi- 
mately 3.5 GPa. Nevertheless, the slope of the calculated 
coexistence line (2.8xl0 6 Pa K^ 1 ) agrees very well with 
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the experimental value (2.7-3.1 xlO 6 Pa K" 1 ) 

At temperatures below 1000 K, the quantum coexis- 
tence curve flattens out and deviates from the straight 
classical line (FIG. [2]). At K, the quantum transi- 
tion pressure is 0.8 GPa higher than the correspond- 
ing classical value. Analysis of our data shows that this 
shift is a direct consequence of the diamond zero-point 
energy being larger than that of graphite. The shape 
of the calculated quantum coexistence line closely re- 
sembles the shape of the Berman-Simon curve obtained 
from experimental thermodynamic properties of diamond 
and graphite [13, [39| . The K transition pressure pre- 
dicted by both the NN and PBE functional (4.7 GPa) 
is again overestimated by 3.3 GPa relative to the exper- 
imental value (1.4 GPa) [U H^|. Based on this obser- 
vation we infer that the positive 3.3 GPa shift of the 
calculated coexistence line is caused by inaccuracies of 
the ab initio PES and not by errors in the NN map- 
ping. This systematic shift is a result of the inability of 
the PBE functional to capture precisely the small differ- 
ence between the energies of graphite and diamond (i.e. 
AU^J 3 ^ = 68 meV/atom is smaller than the average er- 
ror of the PBE functional, 160 meV/atom [I(|). Despite 
this error, the PBE functional and NN predict the zero- 
point energy contributions for diamond and graphite cor- 
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FIG. 2: (color) Graphite-diamond coexistence line. NN re- 
sults are denoted by red, LCBOP+ data by blue, and exper- 
imental data by green color, respectively. 



rectly and, therefore, accurately describe the flattening 
of the coexistence line at the low temperatures. The in- 
set of FIG. [2] shows that the Berman-Simon curve pjjj ]. 
the coexistence line of Bundy 37[ and the experimental 
estimate of the graphite-diamond-liquid triple point 41 1 
are well reproduced in our calculations if the quantum 
NN curve is shifted down by 3.3 GPa to match the ex- 
perimental K transition pressure. 

In the 1000-3000 K range, the coexistence line pre- 
dicted with the LCPOBI+ potential ^ lies -2 GPa 
closer to the experimental line than the NN curve. How- 
ever, the LCBOPI+ potential incorrectly predicts an in- 
crease in the slope of the line below 1000 K and above 
3000 K. As a consequence, the LCBOPI+ triple point 
lies ~4 GPa above the experimental value even though 
the K transition pressure is correctly estimated by 
LCBOPI+. 

In summary, we have demonstrated that despite the 
distinct nature of bonding in graphite and diamond the 
newly developed NN potential predicts numerous proper- 
ties of both phases in quantitative agreement with DFT 
and experimental data. The computational efficiency 
of the NN potential enables an MD study of graphite- 
diamond coexistence of unprecedented accuracy. Com- 
parison of the coexistence lines determined in quantum 
and classical simulations has shown that nuclear quan- 
tum effects are responsible for the experimentally ob- 
served flattening of the coexistence curve at temperatures 
below 1000 K. A detailed MD study of the mechanism of 
the graphite-to-diamond transformation and refinement 
of the NN potential so as to include high-pressure solid 
and liquid phases of carbon are useful follow-on develop- 
ments of this work. 
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